home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / aminet / misc / sci / ephem_src_4_28.lha / objx.c < prev    next >
C/C++ Source or Header  |  1992-05-28  |  30KB  |  1,188 lines
  1. /* functions to save the user-definable objects, referred to as "x" and "y".
  2.  * this way, once defined, the objects can be quieried for position just like
  3.  * the other bodies, with obj_cir(). 
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include <ctype.h>
  9. #ifdef VMS
  10. #include <stdlib.h>
  11. #endif
  12. #include "astro.h"
  13. #include "circum.h"
  14. #include "screen.h"
  15.  
  16. extern char *strcat(), *strcpy(), *strncpy(),
  17. #ifdef AMIGA
  18.     *wgetenv();
  19. #else
  20.     *getenv();
  21. #endif
  22.  
  23. void obj_on(), obj_off(), obj_dfixed(), obj_delliptical(), obj_dhyperbolic(),
  24.     obj_dparabolic(), obj_cir(), crack_year();
  25.  
  26. static char *dbfile;            /* !0 if set by -d option */
  27. #ifdef AMIGA
  28. static char dbfdef[] = "S:ephem.db";     /* default database file name */
  29. #else
  30. static char dbfdef[] = "ephem.db";     /* default database file name */
  31. #endif
  32.  
  33. /* structures to describe objects of various types.
  34.  */
  35. #define    MAXNM        16    /* longest allowed object name, inc \0 */
  36. typedef struct {
  37.     double m_m1, m_m2;    /* either g/k or H/G, depending on... */
  38.     int m_whichm;    /* one of MAG_gk or MAG_HG */
  39. } Mag;
  40. typedef struct {
  41.     double f_ra;    /* ra, rads, at given epoch */
  42.     double f_dec;    /* dec, rads, at given epoch */
  43.     double f_mag;    /* visual magnitude */
  44.     double f_siz;    /* angular size, in arc seconds */
  45.     double f_epoch;    /* the given epoch, as an mjd */
  46.     char   f_name[MAXNM]; /* name */
  47. } ObjF;            /* fixed object */
  48. typedef struct {
  49.     double e_inc;    /* inclination, degrees */
  50.     double e_Om;    /* longitude of ascending node, degrees */
  51.     double e_om;    /* argument of perihelion, degress */
  52.     double e_a;        /* mean distance, aka, semi-maj axis, in AU */
  53.     double e_n;        /* daily motion, degrees/day */
  54.     double e_e;        /* eccentricity */
  55.     double e_M;        /* mean anomaly, ie, degrees from perihelion at... */
  56.     double e_cepoch;    /* epoch date (M reference), as an mjd */
  57.     double e_epoch;    /* equinox year (inc/Om/om reference), as an mjd */
  58.     Mag    e_mag;    /* magnitude */
  59.     double e_siz;    /* angular size, in arc seconds at 1 AU */
  60.     char   e_name[MAXNM]; /* name */
  61. } ObjE;            /* object in heliocentric elliptical orbit */
  62. typedef struct {
  63.     double h_ep;    /* epoch of perihelion, as an mjd */
  64.     double h_inc;    /* inclination, degs */
  65.     double h_Om;    /* longitude of ascending node, degs */
  66.     double h_om;    /* argument of perihelion, degs. */
  67.     double h_e;        /* eccentricity */
  68.     double h_qp;    /* perihelion distance, AU */
  69.     double h_epoch;    /* equinox year (inc/Om/om reference), as an mjd */
  70.     double h_g, h_k;    /* magnitude model coefficients */
  71.     double h_siz;    /* angular size, in arc seconds at 1 AU */
  72.     char   h_name[MAXNM]; /* name */
  73. } ObjH;            /* object in heliocentric parabolic trajectory */
  74. typedef struct {
  75.     double p_ep;    /* epoch of perihelion, as an mjd */
  76.     double p_inc;    /* inclination, degs */
  77.     double p_qp;    /* perihelion distance, AU */
  78.     double p_om;    /* argument of perihelion, degs. */
  79.     double p_Om;    /* longitude of ascending node, degs */
  80.     double p_epoch;    /* reference epoch, as an mjd */
  81.     double p_g, p_k;    /* magnitude model coefficients */
  82.     double p_siz;    /* angular size, in arc seconds at 1 AU */
  83.     char   p_name[MAXNM]; /* name */
  84. } ObjP;            /* object in heliocentric parabolic trajectory */
  85.  
  86. typedef struct {
  87.     int  o_type;    /* current object type; see flags, below */
  88.     int  o_on;        /* !=0 while current object is active */
  89.     ObjF o_f;        /* the fixed object */
  90.     ObjE o_e;        /* the elliptical orbit object */
  91.     ObjH o_h;        /* the hyperbolic orbit object */
  92.     ObjP o_p;        /* the parabolic orbit object */
  93. } Obj;
  94.  
  95. /* o_type */
  96. #define    FIXED        1
  97. #define    ELLIPTICAL    2
  98. #define    HYPERBOLIC    3
  99. #define    PARABOLIC    4
  100.  
  101. /* m_whichm */
  102. #define    MAG_HG        0    /* using 0 makes HG the initial default */
  103. #define    MAG_gk        1
  104.  
  105. static Obj objx;
  106. static Obj objy;
  107.  
  108. #define    DY    0        /* decimal year flag for set_year() */
  109. #define    YMD    1        /* year/mon/day flag for set_year() */
  110.  
  111. /* run when Objx or y is picked from menu.
  112.  * we tell which by the planet code.
  113.  * let op define object and turn it on and off.
  114.  */
  115. void obj_setup(p)
  116. int p;
  117. {
  118.     static char *pr[6] = { /* leave a slot for "On"/"Off" */
  119.         "Fixed", "Elliptical", "Hyperbolic", "Parabolic", "Lookup"
  120.     };
  121.     int f;
  122.     Obj *op;
  123.  
  124.     op = (p == OBJX) ? &objx : &objy;
  125.  
  126.     rechk:
  127.     /* map o_type to popup choice.
  128.      */
  129.     switch (op->o_type) {
  130.     case FIXED: f = 0; break;
  131.     case ELLIPTICAL: f = 1; break;
  132.     case HYPERBOLIC: f = 2; break;
  133.     case PARABOLIC: f = 3; break;
  134.     default: f = 4; break;
  135.     }
  136.  
  137.     ask:
  138.     pr[5] = op->o_on ? "On" : "Off";
  139.     switch (f = popup (pr, f, 6)) {
  140.     case 0: obj_dfixed(op, 0, (char**)0); goto ask;
  141.     case 1: obj_delliptical(op, 0, (char**)0); goto ask;
  142.     case 2: obj_dhyperbolic(op, 0, (char**)0); goto ask;
  143.     case 3: obj_dparabolic(op, 0, (char**)0); goto ask;
  144.     case 4: if (obj_filelookup(p, (char *)0) == 0) obj_on(p); goto rechk;
  145.     case 5: op->o_on ^= 1; break;
  146.     }
  147. }
  148.  
  149. /* turn "on" or "off" but don't forget facts about object the object.
  150.  */
  151. void obj_on (p)
  152. int p;
  153. {
  154.     if (p == OBJX)
  155.         objx.o_on = 1;
  156.     else
  157.         objy.o_on = 1;
  158. }
  159.  
  160. void obj_off (p)
  161. int p;
  162. {
  163.     if (p == OBJX)
  164.         objx.o_on = 0;
  165.     else
  166.         objy.o_on = 0;
  167. }
  168.  
  169. /* return true if object is now on, else 0.
  170.  */
  171. obj_ison(p)
  172. int p;
  173. {
  174.     return ((p == OBJX) ? objx.o_on : objy.o_on);
  175. }
  176.  
  177. /* set an alternate database file name.
  178.  * N.B. we assume the storage pointed to by name is permanent.
  179.  */
  180. void obj_setdbfilename (name)
  181. char *name;
  182. {
  183.     dbfile = name;
  184. }
  185.  
  186. /* fill in info about object x or y.
  187.  * most arguments and conditions are the same as for plans().
  188.  * only difference is that mag is already apparent, not absolute magnitude.
  189.  * this is called by body_cir() for object x and y just like plans() is called
  190.  * for the planets.
  191.  */
  192. void obj_cir (jd, p, lpd0, psi0, rp0, rho0, lam, bet, siz, mag)
  193. double jd;    /* mjd now */
  194. int p;        /* OBJX or OBJY */
  195. double *lpd0;    /* heliocentric longitude, or NOHELIO  */
  196. double *psi0;    /* heliocentric latitude, or 0 if *lpd0 set to NOHELIO */
  197. double *rp0;    /* distance from the sun, or 0 */
  198. double *rho0;    /* true distance from the Earth, or 0 */
  199. double *lam;    /* apparent geocentric ecliptic longitude */
  200. double *bet;    /* apparent geocentric ecliptic latitude */
  201. double *siz;    /* angular size of object, arc seconds */
  202. double *mag;    /* APPARENT magnitude */
  203. {
  204.     Obj *op = (p == OBJX) ? &objx : &objy;
  205.  
  206.     switch (op->o_type) {
  207.     case FIXED: {
  208.         double xr, xd;
  209.         xr = op->o_f.f_ra;
  210.         xd = op->o_f.f_dec;
  211.         if (op->o_f.f_epoch != jd)
  212.         precess (op->o_f.f_epoch, jd, &xr, &xd);
  213.         eq_ecl (jd, xr, xd, bet, lam);
  214.  
  215.         *lpd0 = NOHELIO;
  216.         *psi0 = *rp0 = *rho0 = 0.0;
  217.         *mag = op->o_f.f_mag;
  218.         *siz = op->o_f.f_siz;
  219.         }
  220.         break;
  221.  
  222.     case ELLIPTICAL: {
  223.         /* this is basically the same code as pelement() and plans()
  224.          * combined and simplified for the special case of osculating
  225.          * (unperturbed) elements.
  226.          * inputs have been changed to match the Astronomical Almanac.
  227.          * we have added reduction of elements using reduce_elements().
  228.          */
  229.         double dt, lg, lsn, rsn;
  230.         double nu, ea;
  231.         double ma, rp, lo, slo, clo;
  232.         double inc, psi, spsi, cpsi;
  233.         double y, lpd, rpd, ll, rho, sll, cll;
  234.         double om;        /* arg of perihelion */
  235.         double Om;        /* long of ascending node. */
  236.         double e;
  237.         int pass;
  238.  
  239.         dt = 0;
  240.         sunpos (jd, &lsn, &rsn);
  241.         lg = lsn + PI;
  242.         e = op->o_e.e_e;
  243.  
  244.         for (pass = 0; pass < 2; pass++) {
  245.  
  246.         reduce_elements (op->o_e.e_epoch, jd-dt, degrad(op->o_e.e_inc),
  247.                 degrad (op->o_e.e_om), degrad (op->o_e.e_Om),
  248.                 &inc, &om, &Om);
  249.  
  250.         ma = degrad (op->o_e.e_M
  251.                 + (jd - op->o_e.e_cepoch - dt) * op->o_e.e_n);
  252.         anomaly (ma, e, &nu, &ea);
  253.         rp = op->o_e.e_a * (1-e*e) / (1+e*cos(nu));
  254.         lo = nu + om;
  255.         slo = sin(lo);
  256.         clo = cos(lo);
  257.         spsi = slo*sin(inc);
  258.         y = slo*cos(inc);
  259.         psi = asin(spsi);
  260.         lpd = atan(y/clo)+Om;
  261.         if (clo<0) lpd += PI;
  262.         range (&lpd, 2*PI);
  263.         cpsi = cos(psi);
  264.         rpd = rp*cpsi;
  265.         ll = lpd-lg;
  266.         rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
  267.         dt = rho*5.775518e-3;    /* light travel time, in days */
  268.         if (pass == 0) {
  269.             *lpd0 = lpd;
  270.             *psi0 = psi;
  271.             *rp0 = rp;
  272.             *rho0 = rho;
  273.         }
  274.         }
  275.  
  276.         sll = sin(ll);
  277.         cll = cos(ll);
  278.         if (rpd < rsn)
  279.         *lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
  280.         else
  281.         *lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
  282.         range (lam, 2*PI);
  283.         *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*rsn*sll));
  284.  
  285.         if (op->o_e.e_mag.m_whichm == MAG_HG) {
  286.         /* the H and G parameters from the Astro. Almanac.
  287.          */
  288.         double psi_t, Psi_1, Psi_2, beta;
  289.         beta = acos((rp*rp + rho*rho - rsn*rsn)/ (2*rp*rho));
  290.         psi_t = exp(log(tan(beta/2.0))*0.63);
  291.         Psi_1 = exp(-3.33*psi_t);
  292.         psi_t = exp(log(tan(beta/2.0))*1.22);
  293.         Psi_2 = exp(-1.87*psi_t);
  294.         *mag = op->o_e.e_mag.m_m1 + 5.0*log10(rp*rho)
  295.             - 2.5*log10((1-op->o_e.e_mag.m_m2)*Psi_1
  296.             + op->o_e.e_mag.m_m2*Psi_2);
  297.         } else {
  298.         /* the g/k model of comets */
  299.         *mag = op->o_e.e_mag.m_m1 + 5*log10(rho)
  300.                     + 2.5*op->o_e.e_mag.m_m2*log10(rp);
  301.         }
  302.         *siz = op->o_e.e_siz / rho;
  303.         }
  304.         break;
  305.  
  306.     case HYPERBOLIC: {
  307.         double dt, lg, lsn, rsn;
  308.         double nu, ea;
  309.         double ma, rp, lo, slo, clo;
  310.         double inc, psi, spsi, cpsi;
  311.         double y, lpd, rpd, ll, rho, sll, cll;
  312.         double om;        /* arg of perihelion */
  313.         double Om;        /* long of ascending node. */
  314.         double e;
  315.         double a, n;    /* semi-major axis, mean daily motion */
  316.         int pass;
  317.  
  318.         dt = 0;
  319.         sunpos (jd, &lsn, &rsn);
  320.         lg = lsn + PI;
  321.         e = op->o_h.h_e;
  322.         a = op->o_h.h_qp/(e - 1.0);
  323.         n = .98563/sqrt(a*a*a);
  324.  
  325.         for (pass = 0; pass < 2; pass++) {
  326.  
  327.         reduce_elements (op->o_h.h_epoch, jd-dt, degrad(op->o_h.h_inc),
  328.                 degrad (op->o_h.h_om), degrad (op->o_h.h_Om),
  329.                 &inc, &om, &Om);
  330.  
  331.         ma = degrad ((jd - op->o_h.h_ep - dt) * n);
  332.         anomaly (ma, e, &nu, &ea);
  333.         rp = a * (e*e-1.0) / (1.0+e*cos(nu));
  334.         lo = nu + om;
  335.         slo = sin(lo);
  336.         clo = cos(lo);
  337.         spsi = slo*sin(inc);
  338.         y = slo*cos(inc);
  339.         psi = asin(spsi);
  340.         lpd = atan(y/clo)+Om;
  341.         if (clo<0) lpd += PI;
  342.         range (&lpd, 2*PI);
  343.         cpsi = cos(psi);
  344.         rpd = rp*cpsi;
  345.         ll = lpd-lg;
  346.         rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
  347.         dt = rho*5.775518e-3;    /* light travel time, in days */
  348.         if (pass == 0) {
  349.             *lpd0 = lpd;
  350.             *psi0 = psi;
  351.             *rp0 = rp;
  352.             *rho0 = rho;
  353.         }
  354.         }
  355.  
  356.         sll = sin(ll);
  357.         cll = cos(ll);
  358.         if (rpd < rsn)
  359.         *lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
  360.         else
  361.         *lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
  362.         range (lam, 2*PI);
  363.         *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*rsn*sll));
  364.  
  365.         *mag = op->o_h.h_g + 5*log10(rho) + 2.5*op->o_h.h_k*log10(rp);
  366.         *siz = op->o_h.h_siz / rho;
  367.         }
  368.         break;
  369.  
  370.     case PARABOLIC: {
  371.         double inc, om, Om;
  372.         double lpd, psi, rp, rho;
  373.         double dt;
  374.         int pass;
  375.  
  376.         /* two passes to correct lam and bet for light travel time. */
  377.         dt = 0.0;
  378.         for (pass = 0; pass < 2; pass++) {
  379.         reduce_elements (op->o_p.p_epoch, jd-dt, degrad(op->o_p.p_inc),
  380.             degrad(op->o_p.p_om), degrad(op->o_p.p_Om), &inc, &om, &Om);
  381.         comet (jd-dt, op->o_p.p_ep, inc, om, op->o_p.p_qp, Om,
  382.                     &lpd, &psi, &rp, &rho, lam, bet);
  383.         if (pass == 0) {
  384.             *lpd0 = lpd;
  385.             *psi0 = psi;
  386.             *rp0 = rp;
  387.             *rho0 = rho;
  388.         }
  389.         dt = rho*5.775518e-3;    /* au to light-days */
  390.         }
  391.         *mag = op->o_p.p_g + 5*log10(rho) + 2.5*op->o_p.p_k*log10(rp);
  392.         *siz = op->o_p.p_siz / rho;
  393.         }
  394.         break;
  395.  
  396.     default:
  397.         f_msg ((p == OBJX) ? "Obj X is not defined"
  398.                    : "Obj Y is not defined");
  399.         break;
  400.     }
  401. }
  402.  
  403. /* define obj based on the ephem.db line, s.
  404.  * p is one of OBJX or OBJY.
  405.  * format: name,type,[other fields, as per corresponding ObjX typedef]
  406.  * N.B. we replace all ',' within s with '\0' IN PLACE.
  407.  * return 0 if ok, else print reason why not with f_msg() and return -1.
  408.  */
  409. obj_define (p, s)
  410. int p;    /* OBJX or OBJY */
  411. char *s;
  412. {
  413. #define    MAXARGS    20
  414.     char *av[MAXARGS];    /* point to each field for easy reference */
  415.     char c;
  416.     int ac;
  417.     Obj *op = (p == OBJX) ? &objx : &objy;
  418.  
  419.     /* parse into comma separated fields */
  420.     ac = 0;
  421.     av[0] = s;
  422.     do {
  423.         c = *s++;
  424.         if (c == ',' || c == '\0') {
  425.         s[-1] = '\0';
  426.         av[++ac] = s;
  427.         }
  428.     } while (c);
  429.  
  430.     if (ac < 2) {
  431.         char buf[NC];
  432.         if (ac > 0)
  433.         (void) sprintf (buf, "No type for Object %s", av[0]);
  434.         else
  435.         (void) sprintf (buf, "No fields in %s", s);
  436.         f_msg (buf);
  437.         return (-1);
  438.     }
  439.  
  440.     /* switch out on type of object - the second field */
  441.     switch (av[1][0]) {
  442.     case 'f':
  443.         if (ac != 6 && ac != 7) {
  444.         char buf[NC];
  445.         (void) sprintf(buf,
  446.             "Need ra,dec,mag,D[,siz] for fixed object %s", av[0]);
  447.         f_msg (buf);
  448.         return (-1);
  449.         }
  450.         obj_dfixed (op, ac, av);
  451.         break;
  452.  
  453.     case 'e':
  454.         if (ac != 13 && ac != 14) {
  455.         char buf[NC];
  456.         (void) sprintf (buf,
  457.             "Need i,O,o,a,n,e,M,E,D,H/g,G/k[,siz] for elliptical object %s",
  458.                                     av[0]);
  459.         f_msg (buf);
  460.         return (-1);
  461.         }
  462.         obj_delliptical (op, ac, av);
  463.         break;
  464.  
  465.     case 'h':
  466.         if (ac != 11 && ac != 12) {
  467.         char buf[NC];
  468.         (void) sprintf (buf,
  469.             "Need T,i,O,o,e,q,D,g,k[,siz] for hyperbolic object %s", av[0]);
  470.         f_msg (buf);
  471.         return (-1);
  472.         }
  473.         obj_dhyperbolic (op, ac, av);
  474.         break;
  475.  
  476.     case 'p':
  477.         if (ac != 10 && ac != 11) {
  478.         char buf[NC];
  479.         (void) sprintf (buf,
  480.             "Need T,i,o,q,O,D,g,k[,siz] for parabolic object %s", av[0]);
  481.         f_msg (buf);
  482.         return (-1);
  483.         }
  484.         obj_dparabolic (op, ac, av);
  485.         break;
  486.  
  487.     default: {
  488.         char buf[NC];
  489.         (void) sprintf (buf, "Unknown type for Object %s: %s",
  490.                                 av[0], av[1]);
  491.         f_msg (buf);
  492.         return (-1);
  493.         }
  494.     }
  495.  
  496.     return (0);
  497. }
  498.  
  499. /* if name, then look it up in the ephem database file and set p.
  500.  * else display a table of all objects and let op pick one.
  501.  * p is either OBJX or OBJY.
  502.  * if -d was used use it; else if EPHEMDB env set use it, else use default.
  503.  * return 0 if successfully set object p, else -1.
  504.  */
  505. obj_filelookup (p, name)
  506. int p;            /* OBJX or OBJY */
  507. char *name;
  508. {
  509. /* redefine RCTN,NTR,NTC for column-major order if you prefer */
  510. #define    NLR        (NR-1)        /* number of rows of names.
  511.                      * leave 1 for prompt
  512.                      */
  513. #define    LCW        9        /* screen columns per name */
  514. #define    NLC        9        /* total number of name columns */
  515. #define    NL        (NLR*NLC)    /* total number of names per screen */
  516. #define    RCTN(r,c)    ((r)*NLC+(c))    /* row/col to index */
  517. #define    NTR(n)        ((n)/NLC)    /* index to row */
  518. #define    NTC(n)        ((n)%NLC)    /* index to col (0 based) */
  519.                     /* N.B. all these are 0-based */
  520.     static char prompt[] =
  521.             "RETURN to select, p/n for previous/next page, q to quit";
  522.     FILE *fp;
  523.     char *fn;
  524.     int i, pgn;    /* index on current screen, current page number */
  525.     int r, c;
  526.     char buf[160];    /* longer than any one database line */
  527.     char pb[NC];    /* prompt buffer */
  528.     int readahd;    /* 1 if buffer set from previous loop */
  529.     int choice;    /* index to selection; -1 until set */
  530.     int roaming;    /* 1 while just roaming around screen */
  531.     int abandon;    /* 1 if decide to not pick afterall */
  532.  
  533.     /* open the database file */
  534.     if (dbfile)
  535.         fn = dbfile;
  536.     else {
  537. #ifdef AMIGA
  538.         fn = wgetenv ("EPHEMDB");
  539. #else
  540.         fn = getenv ("EPHEMDB");
  541. #endif
  542.         if (!fn)
  543.         fn = dbfdef;
  544.     }
  545.     fp = fopen (fn, "r");
  546.     if (!fp) {
  547.         (void) sprintf (buf, "Can not open database file %s", fn);
  548.         f_msg(buf);
  549.         return(-1);
  550.     }
  551.  
  552.     /* name is specified so just search for it without any op interaction */
  553.     if (name) {
  554.         int nl = strlen (name);
  555.         int ret = 0;
  556.         while (nxt_db(buf, sizeof(buf), fp) == 0 && strncmp(buf, name, nl))
  557.         continue;
  558.         if (feof(fp)) {
  559.         (void) sprintf (buf, "Object %s not found", name);
  560.         f_msg (buf);
  561.         ret = -1;
  562.         } else
  563.         (void) obj_define (p, buf);
  564.         (void) fclose (fp);
  565.         return (ret);
  566.     }
  567.         
  568.     pgn = 0;
  569.     readahd = 0;
  570.     choice = -1;
  571.     abandon = 0;
  572.  
  573.     /* continue until a choice is made or op abandons the attempt */
  574.     do {
  575.         /* put up next screen full of names.
  576.          * leave top row open for messages.
  577.          */
  578.         c_erase();
  579.         for (i = 0; i < NL; )
  580.         if (readahd || nxt_db (buf, sizeof(buf), fp) == 0) {
  581.             char objname[LCW];
  582.             int ii;
  583.             for (ii = 0; ii < sizeof(objname)-1; ii++)
  584.             if ((objname[ii] = buf[ii]) == ',')
  585.                 break;
  586.             objname[ii] = '\0';
  587.             if (i == NL-1)
  588.             objname[LCW-2] = '\0'; /* avoid scroll in low-r corner*/
  589.             f_string (NTR(i)+2, NTC(i)*LCW+1, objname);
  590.             i++;
  591.             readahd = 0;
  592.         } else
  593.             break;
  594.  
  595.         /* read another to check for eof. if valid, set readahd for next
  596.          * time.
  597.          */
  598.         if (nxt_db (buf, sizeof(buf), fp) == 0)
  599.         readahd = 1;
  600.  
  601.         /* let op pick one. set cursor on first one.
  602.          * remember these r/c are 0-based, but c_pos() is 1-based 
  603.          */
  604.         (void) sprintf (pb, "Page %d%s. %s", pgn+1,
  605.                         feof(fp) ? " (last)" : "", prompt);
  606.         f_prompt(pb);
  607.         r = c = 0;
  608.         roaming = 1;
  609.         do {
  610.         c_pos (r+2, c*LCW+1);
  611.         switch (read_char()) {
  612.         case 'h': /* left */
  613.             if (c == 0) c = NLC;
  614.             c -= 1;
  615.             if (RCTN(r,c) >= i)
  616.             c = NTC(i-1);
  617.             break;
  618.         case 'j': /* down */
  619.             if (++r == NLR) r = 0;
  620.             if (RCTN(r,c) >= i)
  621.             r = 0;
  622.             break;
  623.         case 'k': /* up */
  624.             if (r == 0) r = NLR;
  625.             r -= 1;
  626.             while (RCTN(r,c) >= i)
  627.             r -= 1;
  628.             break;
  629.         case 'l': /* right */
  630.             if (++c == NLC) c = 0;
  631.             if (RCTN(r,c) >= i)
  632.             c = 0;
  633.             break;
  634.         case REDRAW:
  635.             /* start over and skip over prior pages' entries */
  636.             rewind(fp);
  637.             for (i = 0; i < NL*pgn; i++)
  638.             (void) nxt_db (buf, sizeof(buf), fp);
  639.             readahd = 0;
  640.             roaming = 0;
  641.             break;
  642.         case 'p':
  643.             /* if not at first page, start over and skip back one
  644.              * pages' entries
  645.              */
  646.             if (pgn > 0) {
  647.             rewind(fp);
  648.             pgn--;
  649.             for (i = 0; i < NL*pgn; i++)
  650.                 (void) nxt_db (buf, sizeof(buf), fp);
  651.             readahd = 0;
  652.             roaming = 0;
  653.             }
  654.             break;
  655.         case 'n':
  656.             /* if not already at eof, we can go ahead another page */
  657.             if (!feof (fp)) {
  658.             pgn++;
  659.             roaming = 0;
  660.             }
  661.             break;
  662.         case END:
  663.             abandon = 1;
  664.             roaming = 0;
  665.             break;
  666.         case ' ': case '\r':
  667.             choice = NL*pgn + RCTN(r,c);
  668.             roaming = 0;
  669.             break;
  670.         }
  671.         } while (roaming);
  672.     } while (choice < 0 && !abandon);
  673.  
  674.     if (choice >= 0) {
  675.         /* skip first choice entries; selection is the next one */
  676.         (void) rewind (fp);
  677.         for (i = 0; i < choice; i++)
  678.         (void) nxt_db (buf, sizeof(buf), fp);
  679.         (void) nxt_db (buf, sizeof(buf), fp);
  680.         (void) obj_define (p, buf);
  681.     }
  682.     (void) fclose (fp);
  683.     redraw_screen (2);
  684.     return (choice >= 0 ? 0 : -1);
  685. }
  686.  
  687. /* read database file fp and put next valid entry into buf.
  688.  * return 0 if ok, else -1
  689.  */
  690. static
  691. nxt_db (buf, blen, fp)
  692. char buf[];
  693. int blen;
  694. FILE *fp;
  695. {
  696.     char s;
  697.     while (1) {
  698.         if (fgets (buf, blen, fp) == 0)
  699.         return (-1);
  700.         s = buf[0];
  701.         if (isalpha(s))
  702.         return (0);
  703.     }
  704. }
  705.  
  706. /* define a fixed object.
  707.  * args in av, in order, are name, type, ra, dec, magnitude, reference epoch
  708.  *   and optional angular size.
  709.  * if av then it is a list of strings to use for each parameter, else must
  710.  * ask for each (but type). the av option is for cracking the ephem.db line.
  711.  * if asking show current settings and leave unchanged if hit RETURN.
  712.  * END aborts without making any more changes.
  713.  * o_type is set to FIXED.
  714.  * N.B. we don't error check av in any way, not even for length.
  715.  */
  716. static void
  717. obj_dfixed (op, ac, av)
  718. Obj *op;
  719. int ac;
  720. char *av[];
  721. {
  722.     char buf[NC];
  723.     char *bp;
  724.     int sts;
  725.  
  726.     op->o_type = FIXED;
  727.  
  728.     if (set_name (av, op->o_f.f_name) < 0)
  729.         return;
  730.  
  731.     if (av) {
  732.         bp = av[2];
  733.         sts = 1;
  734.     } else {
  735.         static char p[] = "RA (h:m:s): (";
  736.         f_prompt (p);
  737.         f_ra (R_PROMPT, C_PROMPT+sizeof(p)-1, op->o_f.f_ra);
  738. #ifdef AMIGA
  739.         cwrite(") ");
  740. #else
  741.         (void) printf (") ");
  742. #endif
  743.         sts = read_line (buf, 8+1);
  744.         if (sts < 0)
  745.         return;
  746.         bp = buf;
  747.     }
  748.     if (sts > 0) {
  749.         int h, m, s;
  750.         f_dec_sexsign (radhr(op->o_f.f_ra), &h, &m, &s);
  751.         f_sscansex (bp, &h, &m, &s);
  752.         sex_dec (h, m, s, &op->o_f.f_ra);
  753.         op->o_f.f_ra = hrrad(op->o_f.f_ra);
  754.     }
  755.  
  756.     if (av) {
  757.         bp = av[3];
  758.         sts = 1;
  759.     } else {
  760.         static char p[] = "Dec (d:m:s): (";
  761.         f_prompt (p);
  762.         f_gangle (R_PROMPT, C_PROMPT+sizeof(p)-1, op->o_f.f_dec);
  763. #ifdef AMIGA
  764.         cwrite(") ");
  765. #else
  766.         (void) printf (") ");
  767. #endif
  768.         sts = read_line (buf, 9+1);
  769.         if (sts < 0)
  770.         return;
  771.         bp = buf;
  772.     }
  773.     if (sts > 0) {
  774.         int dg, m, s;
  775.         f_dec_sexsign (raddeg(op->o_f.f_dec), &dg, &m, &s);
  776.         f_sscansex (bp, &dg, &m, &s);
  777.         sex_dec (dg, m, s, &op->o_f.f_dec);
  778.         op->o_f.f_dec = degrad(op->o_f.f_dec);
  779.     }
  780.  
  781.     if (set_double (av, 4, "Magnitude: ", &op->o_f.f_mag) < 0)
  782.         return;
  783.  
  784.     if (set_year (av, 5,"Reference epoch (UT Date, m/d.d/y or year.d): ",
  785.                             DY, &op->o_f.f_epoch) < 0)
  786.         return;
  787.  
  788.     if (ac == 7 || !av)
  789.         (void) set_double (av, 6, "Angular Size: ", &op->o_f.f_siz);
  790.     else
  791.         op->o_f.f_siz = 0.0;
  792.  
  793. }
  794.  
  795. /* define an object in an elliptical heliocentric orbit.
  796.  * 13 or 14 args in av, in order, are name, type, inclination, longitude of
  797.  *   ascending node, argument of perihelion, mean distance (aka semi-major
  798.  *   axis), daily motion, eccentricity, mean anomaly (ie, degrees from
  799.  *   perihelion), epoch date (ie, time of the mean anomaly value), equinox year
  800.  *   (ie, time of inc/lon/aop), two magnitude coefficients and optional size.
  801.  * the mag may be H/G or g/k model, set by leading g or H (use H/G if none).
  802.  * if av then it is a list of strings to use for each parameter, else must
  803.  * ask for each. the av option is for cracking the ephem.db line.
  804.  * if asking show current settings and leave unchanged if hit RETURN.
  805.  * END aborts without making any more changes.
  806.  * o_type is set to ELLIPTICAL.
  807.  * N.B. we don't error check av in any way, not even for length.
  808.  */
  809. static void
  810. obj_delliptical(op, ac, av)
  811. Obj *op;
  812. int ac;
  813. char *av[];
  814. {
  815.     op->o_type = ELLIPTICAL;
  816.  
  817.     if (set_name (av, op->o_e.e_name) < 0)
  818.         return;
  819.  
  820.     if (set_double (av, 2, "Inclination (degs):", &op->o_e.e_inc) < 0)
  821.         return;
  822.  
  823.     if (set_double (av, 3, "Longitude of ascending node (degs):",
  824.                 &op->o_e.e_Om) < 0)
  825.         return;
  826.  
  827.     if (set_double (av, 4, "Argument of Perihelion (degs):",
  828.                 &op->o_e.e_om) < 0)
  829.         return;
  830.  
  831.     if (set_double (av, 5, "Mean distance (AU):", &op->o_e.e_a) < 0)
  832.         return;
  833.  
  834.     if (set_double (av, 6, "Daily motion (degs/day):", &op->o_e.e_n) < 0)
  835.         return;
  836.  
  837.     if (set_double (av, 7, "Eccentricity:", &op->o_e.e_e) < 0)
  838.         return;
  839.  
  840.     if (set_double (av, 8, "Mean anomaly (degs):", &op->o_e.e_M) < 0)
  841.         return;
  842.  
  843.     if (set_year (av, 9, "Epoch date (UT Date, m/d.d/y or year.d): ",
  844.                             YMD, &op->o_e.e_cepoch) < 0)
  845.         return;
  846.  
  847.     if (set_year (av, 10, "Equinox year (UT Date, m/d.d/y or year.d): ",
  848.                             DY, &op->o_e.e_epoch) < 0)
  849.         return;
  850.  
  851.     if (set_mag (av, 11, &op->o_e.e_mag) < 0)
  852.         return;
  853.  
  854.     if (ac == 14 || !av)
  855.         (void) set_double (av, 13, "Angular Size @ 1 AU: ", &op->o_e.e_siz);
  856.     else
  857.         op->o_e.e_siz = 0.0;
  858.  
  859. }
  860.  
  861. /* define an object in heliocentric hyperbolic orbit.
  862.  * 11 or 12 args in av, in order, are name, type, epoch of perihelion,
  863.  *   inclination, longitude of ascending node, argument of perihelion,
  864.  *    eccentricity, perihelion distance, reference epoch, absolute magnitude
  865.  *    and luminosity index, and optional size.
  866.  * if av then it is a list of strings to use for each parameter, else must
  867.  * ask for each. the av option is for cracking the ephem.db line.
  868.  * if asking show current settings and leave unchanged if hit RETURN.
  869.  * END aborts without making any more changes.
  870.  * o_type is set to HYPERBOLIC.
  871.  * N.B. we don't error check av in any way, not even for length.
  872.  */
  873. static void
  874. obj_dhyperbolic (op, ac, av)
  875. Obj *op;
  876. int ac;
  877. char *av[];
  878. {
  879.     op->o_type = HYPERBOLIC;
  880.  
  881.     if (set_name (av, op->o_h.h_name) < 0)
  882.         return;
  883.  
  884.     if (set_year(av,2,"Epoch of perihelion (UT Date, m/d.d/y or year.d): ",
  885.                             YMD, &op->o_h.h_ep) < 0)
  886.         return;
  887.  
  888.     if (set_double (av, 3, "Inclination (degs):", &op->o_h.h_inc) < 0)
  889.         return;
  890.  
  891.     if (set_double (av, 4,
  892.         "Longitude of ascending node (degs):", &op->o_h.h_Om) < 0)
  893.         return;
  894.  
  895.     if (set_double(av,5,"Argument of perihelion (degs):", &op->o_h.h_om) <0)
  896.         return;
  897.  
  898.     if (set_double (av, 6, "Eccentricity:", &op->o_h.h_e) < 0)
  899.         return;
  900.  
  901.     if (set_double (av, 7, "Perihelion distance (AU):", &op->o_h.h_qp) < 0)
  902.         return;
  903.  
  904.     if (set_year (av, 8, "Reference epoch (UT Date, m/d.d/y or year.d): ",
  905.                             DY, &op->o_h.h_epoch) < 0)
  906.         return;
  907.  
  908.     if (set_double (av, 9, "g:", &op->o_h.h_g) < 0)
  909.         return;
  910.  
  911.     if (set_double (av, 10, "k:", &op->o_h.h_k) < 0)
  912.         return;
  913.  
  914.     if (ac == 12 || !av)
  915.         (void) set_double (av, 11, "Angular Size @ 1 AU: ", &op->o_h.h_siz);
  916.     else
  917.         op->o_h.h_siz = 0.0;
  918. }
  919.  
  920. /* define an object in heliocentric parabolic orbit.
  921.  * 10 or 11 args in av, in order, are name, type, epoch of perihelion,
  922.  *   inclination, argument of perihelion, perihelion distance, longitude of
  923.  *   ascending node, reference epoch, absolute magnitude and luminosity index,
  924.  *   and optional size.
  925.  * if av then it is a list of strings to use for each parameter, else must
  926.  * ask for each. the av option is for cracking the ephem.db line.
  927.  * if asking show current settings and leave unchanged if hit RETURN.
  928.  * END aborts without making any more changes.
  929.  * o_type is set to PARABOLIC.
  930.  * N.B. we don't error check av in any way, not even for length.
  931.  */
  932. static void
  933. obj_dparabolic(op, ac, av)
  934. Obj *op;
  935. int ac;
  936. char *av[];
  937. {
  938.     op->o_type = PARABOLIC;
  939.  
  940.     if (set_name (av, op->o_p.p_name) < 0)
  941.         return;
  942.  
  943.     if (set_year(av,2,"Epoch of perihelion (UT Date, m/d.d/y or year.d): ",
  944.                             YMD, &op->o_p.p_ep) < 0)
  945.         return;
  946.  
  947.     if (set_double (av, 3, "Inclination (degs):", &op->o_p.p_inc) < 0)
  948.         return;
  949.  
  950.     if (set_double(av,4,"Argument of perihelion (degs):", &op->o_p.p_om) <0)
  951.         return;
  952.  
  953.     if (set_double (av, 5, "Perihelion distance (AU):", &op->o_p.p_qp) < 0)
  954.         return;
  955.  
  956.     if (set_double (av, 6,
  957.         "Longitude of ascending node (degs):", &op->o_p.p_Om) < 0)
  958.         return;
  959.  
  960.     if (set_year (av, 7, "Reference epoch (UT Date, m/d.d/y or year.d): ",
  961.                             DY, &op->o_p.p_epoch) < 0)
  962.         return;
  963.  
  964.     if (set_double (av, 8, "g:", &op->o_p.p_g) < 0)
  965.         return;
  966.  
  967.     if (set_double (av, 9, "k:", &op->o_p.p_k) < 0)
  968.         return;
  969.  
  970.     if (ac == 11 || !av)
  971.         (void) set_double (av, 10, "Angular Size @ 1 AU: ", &op->o_p.p_siz);
  972.     else
  973.         op->o_p.p_siz = 0.0;
  974. }
  975.  
  976. static
  977. set_double (av, vn, pr, fp)
  978. char *av[];    /* arg list */
  979. int vn;        /* which arg */
  980. char *pr;    /* prompt */
  981. double *fp;    /* ptr to double to be set */
  982. {
  983.     int sts;
  984.     char buf[NC];
  985.     char *bp;
  986.  
  987.     if (av) {
  988.         bp = av[vn];
  989.         sts = 1;
  990.     } else {
  991.         f_prompt (pr);
  992.         f_double (R_PROMPT, C_PROMPT+1+strlen(pr), "(%g) ", *fp);
  993.         sts = read_line (buf, 20);
  994.         if (sts < 0)
  995.         return (-1);
  996.         bp = buf;
  997.     }
  998.     if (sts > 0)
  999.         *fp = atof (bp);
  1000.     return (0);
  1001. }
  1002.  
  1003. static
  1004. set_name (av, np)
  1005. char *av[];    /* arg list */
  1006. char *np;    /* name to be set */
  1007. {
  1008.     int sts;
  1009.     char buf[NC];
  1010.     char *bp;
  1011.  
  1012.     if (av) {
  1013.         bp = av[0];
  1014.         sts = 1;
  1015.     } else {
  1016.         (void) sprintf (buf, "Name: (%s) ", np);
  1017.         f_prompt (buf);
  1018.         sts = read_line (buf, MAXNM-1);
  1019.         if (sts < 0)
  1020.         return (-1);
  1021.         bp = buf;
  1022.     }
  1023.     if (sts > 0)
  1024.         (void) strcpy (np, bp);
  1025.     return (0);
  1026. }
  1027.  
  1028. static
  1029. set_year (av, vn, pr, type, yp)
  1030. char *av[];    /* arg list */
  1031. int vn;        /* which arg */
  1032. char *pr;    /* prompt */
  1033. int type;    /* display type: YMD or DY */
  1034. double *yp;    /* ptr to year to be set */
  1035. {
  1036.     int sts;
  1037.     char buf[NC];
  1038.     char *bp;
  1039.  
  1040. #ifdef AMIGA
  1041.     char bb[100];
  1042. #endif
  1043.  
  1044.     if (av) {
  1045.         bp = av[vn];
  1046.         sts = 1;
  1047.     } else {
  1048.         f_prompt (pr);
  1049.         if (type == DY) {
  1050.         double y;
  1051.         mjd_year (*yp, &y);
  1052. #ifdef AMIGA
  1053.         sprintf(bb, "(%g) ", y);
  1054.         cwrite(bb);
  1055. #else
  1056.         (void) printf ("(%g) ", y);
  1057. #endif
  1058.         } else {
  1059.         int m, y;
  1060.         double d;
  1061.         mjd_cal (*yp, &m, &d, &y);
  1062. #ifdef AMIGA
  1063.         sprintf(bb, "(%d/%g/%d) ", m, d, y);
  1064.         cwrite(bb);
  1065. #else
  1066.         (void) printf ("(%d/%g/%d) ", m, d, y);
  1067. #endif
  1068.         }
  1069.         sts = read_line (buf, 20);
  1070.         if (sts < 0)
  1071.         return (-1);
  1072.         bp = buf;
  1073.     }
  1074.     if (sts > 0)
  1075.         crack_year (bp, yp);
  1076.     return (0);
  1077. }
  1078.  
  1079. /* given either a decimal year (xxxx. something) or a calendar (x/x/x)
  1080.  * convert it to an mjd and store it at *p;
  1081.  */
  1082. static void
  1083. crack_year (bp, p)
  1084. char *bp;
  1085. double *p;
  1086. {
  1087.     if (decimal_year(bp)) {
  1088.         double y = atof (bp);
  1089.         year_mjd (y, p);
  1090.     } else {
  1091.         int m, y;
  1092.         double d;
  1093.         mjd_cal (*p, &m, &d, &y);    /* init with current */
  1094.         f_sscandate (bp, &m, &d, &y);
  1095.         cal_mjd (m, d, y, p);
  1096.     }
  1097. }
  1098.  
  1099. /* read next two args from av and load the magnitude members m_m1 and m_m2.
  1100.  * also set m_whichm to default if this is from the .db file, ie, if av!=0.
  1101.  * #,#     -> model is unchanged
  1102.  * g#,[k]# -> g/k
  1103.  * H#,[G]# -> H/G
  1104.  */
  1105. static
  1106. set_mag (av, vn, mp)
  1107. char *av[];    /* arg list */
  1108. int vn;        /* which arg. we use av[vn] and av[vn+1] */
  1109. Mag *mp;
  1110. {
  1111.     int sts;
  1112.     char buf[NC];
  1113.     char *bp;
  1114.  
  1115.     if (av) {
  1116.         mp->m_whichm = MAG_HG;    /* always the default for the db file */
  1117.         bp = av[vn];
  1118.         sts = 1;
  1119.     } else {
  1120.         /* show both the value and the type of the first mag param,
  1121.          * as well as a hint as to how to set the type if desired.
  1122.          */
  1123.         (void) sprintf (buf, "%c: (%g) (g# H# or #) ",
  1124.                 mp->m_whichm == MAG_HG ? 'H' : 'g', mp->m_m1);
  1125.         f_prompt (buf);
  1126.         sts = read_line (buf, 9);
  1127.         if (sts < 0)
  1128.         return (-1);
  1129.         bp = buf;
  1130.     }
  1131.     if (sts > 0) {
  1132.         switch (bp[0]) {
  1133.         case 'g':
  1134.         mp->m_whichm = MAG_gk;
  1135.         bp++;
  1136.         break;
  1137.         case 'H':
  1138.         mp->m_whichm = MAG_HG;
  1139.         bp++;
  1140.         default:
  1141.         /* leave type unchanged if no prefix */
  1142.         break;
  1143.         }
  1144.         mp->m_m1 = atof (bp);
  1145.     }
  1146.  
  1147.     if (av) {
  1148.         bp = av[vn+1];
  1149.         sts = 1;
  1150.     } else {
  1151.         /* can't change the type in the second param */
  1152.         (void) sprintf (buf, "%c: (%g) ",
  1153.                 mp->m_whichm == MAG_HG ? 'G' : 'k', mp->m_m2);
  1154.         f_prompt (buf);
  1155.         sts = read_line (buf, 9);
  1156.         if (sts < 0)
  1157.         return (-1);
  1158.         bp = buf;
  1159.     }
  1160.     if (sts > 0) {
  1161.         int ok = 0;
  1162.         switch (bp[0]) {
  1163.         case 'k':
  1164.         if (mp->m_whichm == MAG_gk) {
  1165.             bp++;
  1166.             ok = 1;
  1167.         }
  1168.         break;
  1169.         case 'G':
  1170.         if (mp->m_whichm == MAG_HG) {
  1171.             bp++;
  1172.             ok = 1;
  1173.         }
  1174.         break;
  1175.         default:
  1176.         ok = 1;
  1177.         break;
  1178.         }
  1179.         if (ok)
  1180.         mp->m_m2 = atof (bp);
  1181.         else {
  1182.         f_msg ("Can't switch magnitude models at second parameter.");
  1183.         return (-1);
  1184.         }
  1185.     }
  1186.     return (0);
  1187. }
  1188.